Dynamical Screening in Correlated Electron Materials 
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We present an efficient method for incorporating the dynamical effects of the screening of the 
Hubbard U by electronic degrees of freedom in the solid into the single site dynamical mean field 
approximation. The formalism is illustrated by model system calculations which capture the es- 
sential features of the frequency dependent interactions proposed for Gd, Ni, SrVOa and other 
compounds. Screening leads to shifts in the metal-insulator phase boundary, changes in the spectral 
function near the Mott-Hubbard gap edge and to a renormalization of the quasiparticle weight. 
Hubbard bands are generically neither separated by the screened nor the unscreened interaction en- 
ergy, implying that the common practice of extracting the Hubbard U from the energies of features 
in photoemission and inverse photoemission spectra requires reexamination. 
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'Strongly correlated electron systems' are a central 
topic in electronic condensed matter physics [l[ . The low 
energy physics of these systems is typically described by 
an effective Hamiltonian which models the behavior of 
relatively localized d or / orbitals and is obtained (at 
least notionally) via a "downfolding" procedure in which 
other degrees of freedom are integrated out. A crucial as- 
pect of the effective Hamiltonian is an interaction which 
acts to suppress local number fluctuations. This inter- 
action is typically parametrized by a number, the "Hub- 
bard U" . However, U is generically dynamical: a density 
fluctuation in a correlated orbital produces electric fields, 
which other degrees of freedom will act to screen, result- 
ing in a frequency-dependent renormalization. Screen- 
ing has been observed as a shift in excitation energies 
in experiments comparing the surface of solid Ceo to Cgo 
films on silver Q and has been computed using variations 
of the 'random phase approximation' (RPA) [3-5]. The 
renormalizations are found to be strong in many cases. 
In Gd, U(uj) rises from a static value of about 6.5 eV 
to about 17 eV as the frequency lu is raised from to 
~ 3 eV, while in Ce, U(ui) changes from ~ 3.5 eV to 
~ 7 eV as u> is increased from to ~ 4 eV 

While the study of the strong correlation effects in- 
duced by an instantaneous interaction is well advanced 
thanks to the development of dynamical mean field the- 
ory (DMFT) our ability to treat frequency depen- 
dent interactions has been limited. The most widely 
used method for solving the DMFT equations has been 
the Hirsch-Fye algorithm Q, which is based on a time- 
discretization and decoupling of interaction terms by 
auxiliary fields. Frequency dependent interactions lead 
to a proliferation of decoupling fields which become pro- 
hibitively expensive to sample. A treatment of screening 
effects within exact diagonalization, numerical renormal- 
ization group and other Hamiltonian based methods re- 
quires the explicit introduction of many bosonic modes, 
leading to a Hilbert space which is too large to be handled 



numerically. Considerations of this sort have led to the 
belief (see e.g. Ref. 0]) that DMFT simulations with fre- 
quency dependent interactions are exceedingly difficult. 

Here we show that this is not the case: the recently 
developed 'hybridization expansion' diagrammatic quan- 
tum Monte Carlo method [1, Q can be used to treat 
models with an arbitrary frequency dependence of the 
on-site repulsion U (w) at negligible additional computa- 
tional cost, opening the door to a systematic investiga- 
tion of screening effects in correlated electron materials. 
We begin our discussion by recalling that the downfoldcd 
models used to describe the correlated degrees of freedom 
in a transition metal or actinide involve a one electron 
part and an interaction part. The parameters describ- 
ing the one-electron physics are obtained by projecting a 
band theory calculation onto a set of distinguished ( "d" ) 
orbitals and are in principle energy dependent. If the 
d orbitals are correctly chosen the energy dependence is 
negligible [13, EH so the one electron part may be mod- 
elled as a tight-binding-likc Hamiltonian -ffband- 

The interaction part is obtained (sec, e.g. Ref. HI) 
by screening the bare Coulomb interaction e 2 /|r — r'\ 
with real and virtual transitions involving the orbitals 
which are integrated out, projecting the result onto the 
d manifold and retaining only the fully site-diagonal 
terms. One finds two kinds of terms: an instantaneous 
interaction and a screening contribution. The instanta- 
neous interaction may be represented as a Hamiltonian 
term H- mt which takes the usual Slatcr-Kanamori form 
Hi n t = J2i \UNf + .... with Ni the number operator for 
electrons in the d manifold on site i and the ellipsis denot- 
ing exchange, 'pair hopping' and other terms which in- 
volve operators such as the spin and angular momentum 
which commute with N. The screening contribution cou- 
ples the site densities at different times and is expressed 
as a contribution Sw = i j drdr' N(t)W(t - t')N(t') 
to the effective action. The screening function W(t) is 
negative, depends only on one time (or frequency) argu- 
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ment and couples only to the site density Ni, reflecting 
the physics of dynamical screening. It is generated by a 
screening spectral function lmW(ui): 
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with W Uo (t) = cosh [(t - § )w ]/sinh [^] for < t < 
(3 and W Uo (t) = Wu (T + nj3) (n integer). Other interac- 
tion parameters acquire only a very weak frequency de- 
pendence Q because they involve operators correspond- 
ing to higher multipoles which are weakly screened. 

We recast the time dependent interaction in Hamilto- 
nian form by using Hubbard-Stratonovich techniques to 
introduce boson operators whose spectrum reproduces 
ImlF. We define A 2 o = — TrImW (uj ) and at each 
site i and each frequency ujq we employ the identity 
("*" denotes integration over time) e 2 A »o Ar,H '"o ,jV ' )l »o = 

$ V^ e-^* W ^*^-^°* N ^° . Noting that = 
(— <9 2 + w 2 )/(2wo) (the periodicity under r — > t + n(3 
implies a derivative discontinuity at r = n/3) and identi- 
fying 9 r Uo as wo times the momentum H Uo conjugate to 
4> UQ we see that the time dependent interaction is equiv- 
alent to a Hamiltonian model in which the site density is 
coupled via A Wo to a set of oscillators, 
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To solve the Hamiltonian H 
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we employ dynamical mean field theory [y, which 
computes the solution to the correlated electron prob- 
lem from the solution of an auxiliary quantum impurity 
model with interactions given by the local interactions 
of the original model and an impurity-bath hybridiza- 
tion term i?hyb = cWd + H.c. with V determined by 
a self-consistency condition. Following Ref. [§| we solve 
the impurity model by expanding the partition function 
in powers V, collecting diagrams into determinants of 
matrices of hybridization functions and using a Monte 
Carlo procedure to sample the resulting sum of determi- 
nants. As in Ref. 12| we treat the electron-boson cou- 
pling by a canonical transformation which shifts (f> Uo by 
\/2 X Uo N/uio , changes the instantaneous interaction U to 
the screened value 



U scr = U + 2 
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and shifts the chemical potential (i to jLt scr = /x + 
Io° Im ^' ,Jo - ) . The transformation also multiplies d, 
rft by e^XT^on^VSA^/^o ( s = 1 for rf t an d s = -1 
for d). The result is that a term in the hybridiza- 
tion expansion with 2n hybridization events at times 
< T\ < T2 < . . . < Tm < /3 acquires an extra weight 
factor w S creen({Ti}) given by an exponential of correlators 
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FIG. 1: Illustration of an order n = 4 diagram for the one 
orbital Hubbard model. Empty (full) circles and squares rep- 
resent V* (V) hybridization events. Dashed lines indicate in- 
teractions K(t) connecting all pairs of hybridization events. 



of nonintcracting boson operators, which evaluates to: 
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It follows from Eqs. Q and ([3]) and the illustration in 
Fig. [1] that the screening contribution amounts to a non- 
local interaction between all pairs of hybridization events. 
Since the computational bottleneck is the handling of the 
hybridization matrix determinants, the additional weight 
factor lOgcreen does not significantly slow down the simu- 
lation, which remains very efficient fl3| . 

We now illustrate the method with single-site, single- 
orbital DMFT calculations for a semi-circular density of 
states of bandwidth D = 4 and inverse temperature 
j3 = 50. All calculations are performed in the para- 
magnetic phase and for half filling. We consider two 
screening functions: (i) K(t) = — (A/wo) 2 (cosh((/3/2 — 
t)u>o) — cosh(/3wo/2))/ sinh(/?wo/2) corresponding to a 
delta- function TmW(cj) = — X 2 tt (<5(w — ujq) — <5(w + wo)) 
and (ii) K{t) = aln[l + /3w c sin(7TT/,£i)/7r)] corresponding 
to "Ohmic" screening ImlF^w) ~ — uttlo at small w, with 
an ultraviolet cutoff at w c . Model (i) implies that the 
screened interaction is U SCI = U— 2X 2 /ujq and RcVF(w) = 
2A 2 wo/(w 2 — Wo); for model (ii) U scr = U — 2aw c and 
KeW(uj) = aw In |(w c + w)/(w c — w)| — 2aw c . In princi- 
ple the coupling strengths can be made large enough to 
drive U SCI negative (overscreening) , but in the absence of 
phonons overscreening is not believed to occur in real ma- 
terials, limiting the range of physically relevant couplings 
to A < y/ujQU/2 or a < U/(2u) c ). Model (i) is a rough 
rcprcsention of results obtained from constrained RPA 
calculations for Gd and SrV03 and model (ii) captures a 
characteristic feature of paramagentic Ni [3|-|5|. 

The calculated metal-insulator phase diagrams for our 
two models are shown in the left-hand panels of Fig [2] 
At our simulation temperature the critical U for the 
metal insulator transition of the unscreened model is 
U C 2{P = 50) w 5.1 and the transition is first-order [f|. 
As expected on physical grounds, increasing the strength 
of the screening shifts the metal-insulator transition to 
larger values of the bare interaction. The transition re- 
mains first order but the coexistence region becomes nar- 
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FIG. 2: Left panels: metal-insulator phase diagram in the 
space of bare interaction U and screening coupling strength A 
(a) for indicated values of the screening frequency uio (w c ). 
Upper panel: model (i); lower panel, model (ii). Dashed 
lines: screening strength at which the screened interaction 
User changes sign. Right panels: phase diagram in the space 
of screened and bare interaction. 



X=2.5, U scr =5.83 



X=3.5, lC=1.83 





rower as the screening effect increases (the two metal- 
Mott insulator phase boundaries are shown for model (i) 
at (jJq = 5; for the other cases we show the stability re- 
gion of the metallic phase). If the coupling strength is 
increased into the unphysical overscreening regime, one 
finds a first order transition to a bipolaronic insulator. 

The right hand panels of Fig. [2] present the phase dia- 
grams as a function of screened interaction and coupling 
strength (parametrized as ratio of unscreened to screened 
interaction). For U /U scr ~ 2-3, typical of values found in 
RPA-type calculations, the critical screened interaction 
is \ to | of the U C 2 defined in the unscreened model. 
The dependence on screening frequency is weak. 

We next consider the effect of screening on the elec- 
tron spectral functions. In the model without screening 
the spectral function is characterized by peaks at uj = 
(if the model is in the metallic phase) and w ~ ±[//2. 
Figure [3] shows that the situation is quite different in the 
screened case. The top panel presents the changes that 
occur as the screening strength is varied at fixed bare 
interaction U — 10 and screening frequency w = 3. A 
multipeak structure is evident, with a broad high energy 
feature at an energy set by a combination of the bare U 
and A, a lower energy sideband and, in the metallic case, 
a (split) peak at ui = 0. While the peak position of the 
lower energy sidebands varies roughly in parallel with the 
screened U its energy is a parameter-dependent fraction 
of U SCT /2 (approximately 60%). 

The lower panel of Fig. [3] presents the variation of the 



FIG. 3: Evolution of the spectral function across the metal- 
insulator transition for model (i). Top panel: effect of the 
variation of the screening strength at fixed U = 10, cjo = 3. 
Bottom panel: effect of the variation of the screening fre- 
quency at fixed U = 8 and U S ct — 3. Light (heavy) dashed 
line: spectrum of the unscreened model for U = 8 (U — 3). In- 
set: frequency dependent interaction KeU(to) = U + KeW (to) . 



spectral function with screening frequency for fixed bare 
and screened interactions U = 8 > U C 2, User — 3 < U C 2- 
For reference the spectrum in the absence of screening 
is also shown. At the lowest value of the screening fre- 
quency, ljq = 1, the model is in its insulating phase and 
the spectrum is a broadened version of that of the un- 
screened insulator. For all other screening frequencies, 
the model is in the metallic phase. We see that the po- 
sitions of both the lower energy sideband and the high 
energy peak increase as the screening frequency is in- 
creased. At the highest screening frequency the lower 
energy portion of the spectrum begins to approach that 
expected in the antiadiabatic limit, with a peak near 
U SCI /2. We therefore interpret the lower energy feature 
as the "screened U" sideband, but note that in general its 
energy does not yield a good measure of U scr . Even at the 
highest screening frequency (u)q — AD — 16), substantial 
spectral weight exists at high energies (w ~ ujq). 

FigureUshows spectra computed for the Ohmic screen- 
ing model. The spectra are qualitatively similar to those 
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FIG. 4: Model (ii): evolution of the spectral function for 
U — 8, t/scr = 3 and indicated values of cj c . Light (heavy) 
dashed line: spectrum of the unscreened model for (7 = 8 
(U = 3). Inset: ReU(u) = U + ReW(u)). 



computed for the plasmon model, but the high energy 
feature appears as a broad tail and not as a separate 
peak. 

In summary, we have described a simple and effi- 
cient algorithm to treat fermionic lattice models with 
arbitrary frequency dependent interactions within single- 
site DMFT. The only restrictions are that the external 
screening degrees of freedom may be represented as non- 
interacting bosons, have an analytically known commu- 
tator with electron creation operators and couple to a 
quantity which is conserved by the on-site Hamiltonian 
(this last restriction prevents a direct application of our 
formalism to clusters). While we have presented results 
for the one-band Hubbard model, we emphasize that our 
method is applicable to multi-orbital models with gen- 
eral Slater-Kamanori interactions. The method there- 
fore opens the door for efficient DMFT simulations of 
the properties of strongly correlated compounds with ar- 
bitrary energy dependence of the interaction parameters, 
as well as self-consistent GW+DMFT calculations. We 
also note that the extended-DMFT variants of dynam- 
ical mean field theory lead to bosonic problems with a 
structure very similar to the problem we have considered 

a 

The frequency dependent U (u>) may be characterized 
by three numbers: a bare (unscreened) interaction, a 
screened interaction, and a screening frequency. If the 
screening frequency is very high compared to both the 
conduction electron bandwidth and the screened inter- 
action, then the physics is well described by an effec- 
tive model with an instantaneous interaction equal to 
the screened i7 scr . If the screening frequency is very low 
compared to these scales, then one recovers the familiar 
electron-phonon physics. We have shown here that the 
crossover between the two regimes is very broad. Over 
wide parameter ranges the dynamical nature of the in- 



teraction is important: a model with effective interaction 
equal to the screened one provides poor estimates of the 
location of the metal-insulator transition, the value of the 
Mott Hubbard gaps and the locations and line shapes of 
the shakeoff features in the spectral function. In partic- 
ular, the peak positions in the spectal functions do not 
provide quantitative estimates of either the screened or 
unscreened U values. 

Our work suggests several directions for future re- 
search. Application of the method to real materials is in 
progress. Our approach may also be useful as a solver for 
E-DMFT calculations. Finally, we note that Assaad and 
Lang have presented a method for treating bosons within 
the weak coupling formalism [lij . While the weak cou- 
pling approach is less efficient than the hybridization ex- 
pansion for single-site (multiorbital) calculations, it can 
be applied to cluster dynamical mean field theory. An 
extension of Assaad's method to investigate the effects 
of screening in clusters would be of interest. 
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